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Abstract. We study tilted perfect fluid cosmological models with a constant equation of state parameter 
i in spatially homogeneous models of Bianchi type VI;, using dynamical systems methods and numerical 

i experimentation, with an emphasis on their future asymptotic evolution. We determine all of the equilibrium 

f - ^ ■ points of the type VI^ state space (which correspond to exact self-similar solutions of the Einstein equations, 

^v^j ' some of which are new), and their stability is investigated. We find that there are vacuum plane- wave 

solutions that act as future attractors. In the parameter space, a 'loophole' is shown to exist in which there 
are no stable equilibrium points. We then show that a Hopf-bifurcation can occur resulting in a stable 
closed orbit (which we refer to as the Mussel attractor) corresponding to points both inside the loophole 
and points just outside the loophole; in the former case the closed curves act as late-time attractors while 
in the latter case these attracting curves will co-exist with attracting equilibrium points. In the special 
Bianchi type III case, centre manifold theory is required to determine the future attractors. Comprehensive 
numerical experiments are carried out to complement and confirm the analytical results presented. We note 
that the Bianchi type VI?, case is of particular interest in that it contains many different subcases which 
exhibit many of the different possible future asymptotic behaviours of Bianchi cosmological models. 
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O ' 1. Introduction 

t> 

In this paper, we shall study tilted perfect fluid cosmological models with a constant equation of state 
' parameter, 7, in spatially homogeneous Bianchi models of type VI^ using the formalism introduced in 

O" 1 DP- We introduce expansion- normalised variables [2], determine the equilibrium points and their stability 
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properties and consequently investigate the future asymptotic behaviour of the models and determine the 
late-time asymptotic states, using dynamical systems methods and numerical experimentation. In spatially 
homogeneous cosmological models the universe is foliated into space-like hypersurfaces (defined by the group 
orbits of the respective model) [5], [31 HI [SI IS] • For these perfect fluid models there are two naturally defined 
time-like vector fields (i.e., congruences): the unit vector field, normal to the group orbits and hence 
orthogonal to the surfaces of transitivity (the 'geometric' congruence), and the four- velocity, w M , of the 
fluid (the 'matter' congruence). If u M is not aligned with ri M , the model is called tilted (and non-tilted or 
orthogonal otherwise) [7]. The geometric congruence is necessarily geodesic, vorticity-free and acceleration- 
free. The matter congruence, on the other hand, is not necessarily geodesic and can have both vorticity and 
acceleration. 

We shall follow convention and use the kinematical quantities associated with the normal congruence 
n M (rather than the fluid flow u M ) of the spatial symmetry surfaces as variables. This avoids the possible 
singular behaviour experienced by the fluid observers in these models El [101 E] (these papers also 
give the appropriate boost-formluae relating the kinematic variables and curvatures in the two frames). In 
principle, it is also necessary to investigate the behaviours in the fluid frame as the asymptotes may differ 
(which is the case for so-called 'whimper singularities' [2]); however, this will not be pursued here (criteria 
for when this may happen are given in We will also focus on the dynamical behaviour rather than the 

physical effects of these models; observational consequences of (tilted) Bianchi models have been studied in, 
for example, (EH O [15] . 

We will assume a perfect-fluid matter source with p = (7 — l)/x as equation of state, where fi is the energy 
density, p is the pressure, and 7 is a constant. Causality then requires 7 to be in the interval < 7 < 2. 
A positive cosmological constant may also be included in the models. Tilted SH cosmologies with a 7-law 
perfect fluid source have been studied by a number of authors. In 1G the stability of non-tilted universes 
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4n these papers the physical properties of these models, and particularly the observations and singularity structure in 
models in which the tilt becomes extreme asymptotically (i.e., v 2 — » 1) as measured by observers moving with the geometric 
congruence, was studied. 
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Table 1 . Previous studies of the dynamical behaviour of tilted Bianchi models 
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General 


32j 



with respect to tilt was studied and Table[T]shows previous studies of tilted universes of Bianchi type II- VIII. 
It has been proven that if the matter obeys the strong energy condition, the positive pressure criteria, the 
dominant energy condition, and a matter regularity condition all Bianchi type IX models recollapse to the 
future QI1[]J3[IS|20]. 

Of the most general ever-expanding Bianchi models (namely VTfc, VII/, and VIII) only the dynamics of 
the type Vlh models remain to be studied. This paper aims to fill this gap in the study of the behaviour 
of general spatially homogeneous Bianchi models. The behaviour of these models has been shown to be 
both interesting and surprising. In particular, for the Bianchi type VIII models (and Bianchi models of type 
VII [31]), the state space is unbounded and consequently, for all non-inflationary perfect fluids, one of the 
curvature variables grows without bound at late times. It was found that in Bianchi type VIII models [32] 
with fluids stiffer than dust (1 < 7 < 2), the fluid will in general tend towards a state of extreme tilt, while 
for fluids less stiff than dust (0 < 7 < 1) the fluid will in the future be asymptotically non-tilted. Using both 
dynamical systems theory and a detailed numerical analysis, the late-time behaviour of tilting perfect fluid 
Bianchi models of types IV and VII?,, was studied, and it was found that the plane waves are the only future 
attracting equilibrium points for non- inflationary fluids in Bianchi type VII,, models 22 . A tiny region of 
parameter space (the "loophole") in the Bianchi type IV model was shown to contain a closed orbit which 
acts as an attractor pQ. From a detailed numerical analysis it was found that at late times the normalised 
energy-density tends to zero and the normalised variables 'freeze' into their asymptotic values, and it was 
then shown that there is an open set of parameter space in the type VII/, models in which solution curves 
approach a compact surface that is topologically a torus [22] . 

In this work we determine all of the equilibrium points of the type VI/, state space. These equilibrium 
points correspond to exact self-similar solutions of the Einstein equations and play a special role in the 
general evolution of the system [33j . In particular, the stability of these solutions are determined. Some of 
these solutions are new (and one is given implicitly). A complete catalogue of self-similar solutions is given 
and the late time attractors are classified. A detailed numerical analysis is included, complementing the 
analytical results. We note that the Bianchi type VI?, case is very complicated, with many cases to consider, 
and the resulting dynamics include many of the different types of behaviour described in previous work. 

This paper is organised as follows. In the following section we present the equations of motion for the 
tilted Bianchi type VI/, models and discuss the state space and various invariant sets of physical interest. In 
section [3] we present monotonic functions and give a list of all equilibrium points. Section [4] is devoted to 
the analysis of the late-time behaviour and in section [5] we give a synopsis of the numerical analysis done. 
The final section is devoted to a discussion of the ramifications of our results. 



2. Equations of motion 
2.1. The orthonormal frame approach. The line-element of a Bianchi cosmology can be written 

(2.1) ds 2 = -dt 2 + 6 ab u a uj b , 
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where t is the co-moving cosmological time. The one-forms us a are left-invariant one-forms on the hypersur- 
faces spanned by the group orbits. They obey the commutator relation 

(2.2) {&u, a )^ = -^C a bc u, b A« c , C a bc = e bcd n da + a b 6 a c - a c S a bl 

where n ab is a symmetric tensor, a c — — (l/2)C a ac and _L means projection onto the spatial hypersurfaces. 
Furthermore, the Jacobi identity implies n bc a c = 0. 

Since n is a symmetric spatial tensor we can characterise it in terms of its eigenvalues, (ni,n 2 , n%). We 
are interested in the Bianchi type Vlh models for which a c ^ 0, and hence, one of the eigenvalues of n ab is 
necessarily zero: n\ — (say). Moreover, for the Bianchi type VI ^ models we have the relation a c a c = ft.TT.2n3, 
which defines the group parameter h. 

The geometric (or normal) congruence, rz/ 1 , is given by n = d/dt. It is also useful to define the shear and 
the Hubble scalar associated with the congruence n^: 

(2.3) H = g^/i) = n w - Hhn V , 

where hp V is the spatial metric on the hypersurfaces spanned by the group orbits. The matter variables are 
chosen to be the energy density, fi, and the tilt-velocity, v a , which is defined as the 3- velocity of the fluid 
with respect to the geometric (or normal) congruence, n^. The equations of motion can now be written 
down in terms of the Hubble scalar, H; the shear, <j ab \ the curvature variables n ab and a c ; and the matter 
variables fi and v a . 

In the dynamical systems approach it is common to introduce expansion-normalised variables (we divide 
the variables with the appropriate powers of H). The paper [T] contains all the details regarding the 
determination of the evolution equations for the tilted cosmological models under consideration. These 
equations, written in gauge invariant form, allow one to choose the gauge that is best suited to the application 
at hand. Here, we shall adopt the ./V-gauge in which the function N x is purely imaginary; this is ensured 
by the choice <fi' = -\/3AE_, where A is defined by N = AIm(N x ). The evolution equation for N can then 
be replaced with an evolution equation for A, which ensures a closed system of equations. For a qualitative 
analysis of these models, the iV-gauge is preferable since the resulting dynamical system is well defined in 
the Bianchi Vlft case, in particular, |A| < 1. The relation between these variables (and their corresponding 
differential equations) and the models they represent are discussed in more detail in [2J. 

In the notation of [T] , the expansion-normalised anisotropy and curvature variables used in this paper are: 

S 1 = S 12 +iSi3, S x =S_+iS 23 , N x =iN. 

We will also adopt the dimensionless time parameter r, which is related to the cosmological time t via 
dt/dr = (1/H), where H is the Hubble scalar. 

Using expansion-normalised variables, the equations of motion in the iV-gauge are (see [I] for the complete 
derivation of the equations) : 

(2.4) S' + 

(2.5) E'_ 

(2.6) Si 2 

(2.7) E' 13 

(2.8) E 23 

(2.9) N' 

(2.10) A' 

(2.11) A' 



(q - 2)E + + 3(E? 2 + E? 3 ) - 27V 2 + JL- {-2v\ + v\ + v 2 3 ) 

(q-2- 2V3S 23 A)S_ + \/3(S 2 2 - E 2 3 ) + 2AN + {v 2 2 - v 2 3 ) 

(q - 2- 3E+ - V3£_) S 12 - V3 (S 23 + S_A) S X3 + ^J^i«2 

(q - 2 - 3S+ + V3S_) E 13 - (S 23 - £-A) S 12 + ^^v lV3 

(q - 2)S 23 - 2V3N 2 \ + 2V3AE2_ + 2V3Si 2 Si 3 + ^P^-v 2 v 3 

G+ 

(q + 2S+ + 2V3S 23 A^ N 

2^E 23 (1 - A 2 ) 
(« + 2E+U. 
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The equations for the fluids are 

n 



(2.12) 


n' 


(2.13) 


v[ 


(2.14) 


v' 2 


(2.15) 


v' 3 


(2.16) 


V 


where 





a fa - (37 - 2) + 2 1 Av 1 + [2,7(7 - 1) - (2 - 7) - ~fS\ V 2 } 
(T + 2E+) vi - 2v / 3Si 3 «3 - 2V3T, 12 v 2 - A («f + «|) - V3A - 



T-£_, 



T - £+ + V3E 

7(1 - 7 2 ) 
1 - (7- 1)V 2 



V3£_ ) v 2 - V3 (£ 23 + E-A) u 3 + \/3AiVwiU3 
V3 — (S23 — £_ A) w 2 — v / 3AA r wiw 2 
[(3 7 - 4) - 2( 7 - l)A Vl - S] , 



A + VSNj 
A - VSN^j 



viv 2 



viv 3 



q 

s 2 = 

s = 

V 2 = 

T = 

G + = 



= 2£ 2 



1 (3 7 - 2) + (2 - j)V 2 



2 

£ 2 _ 



i + ( 7 -i)y 2 

£ 2 



E 2 



^afeC C , C C a 



J 13 
1, «' 



y2 
^23 



'3' 



[(3 7 - 4) - 2(7 - l)Av x ] (1 - I/ 2 ) + (2 - 7 )7 2 <S 



1 - (7- 1)V 2 



1 + (7 - i)v 2 



These variables are subject to the constraints 

(2.17) 1 

(2.18) 



£ 2 + A 2 + N 2 4 
2£+yl + 2£_A 



G 



+ 



(2.19) 

(2.20) 
(2.21) 



= - 



= A 



E 12 (iV + v / 3A) + £i 3 AiV 

Si 3 (JV- V3A) + E 12 XN 
3/i(l-A 2 ) A 2 . 



7^1f 2 
G+ 

7^3 



The parameter 7 will be assumed to be in the interval 7 G (0,2). The generalized Friedmann equa- 
tion (|2.17p yields an expression which effectively defines the energy density fl We will assume that 
this energy density is non-negative: O > 0. Therefore, the state vector can thus be considered X = 
[£ + , £12, £13, £ 23 , A, A, A, ux, «2, W3] modulo the constraint equations (|2.18|) - (|2.21[) . Thus the dimen- 
sion of the physical state space is seven (for a given value of the parameter h). Additional details are 
presented in pQ. 

The dynamical system is invariant under the following discrete symmetries : 

[£+, E_,Ex2, Ei 3 , E 23 , A, A, A, Vi,v 2 ,v 3 ] h-> [E + , E_, E 12 , Ei 3 , E 23 , -A, A, -A, -vi, -v 2 , -v 3 ] 
[E + , £_, £ 12 , Ei 3 , E 23 , A, A, A, v x , v 2 , v 3 ] h-> [E + , -E_, E 13 , E 12 , E 23 , -A", A, A, Vi,v 3 , v 2 ] 
[E+, E_, E12, E 13 , E 23 , A, A, A, vi,v 2 ,v 3 ] h-> [£+,£_, ±£ 12 , =f£ 13 , -E 23 , A, -A, A, v\, ±v 2 ,=pv 3 ] 
[E+, E_, Ex2, Ex3, E 23 , A, A, A, vi,v 2 ,v 3 ] ^ [£+, £_, -Ex2, -£13, £23, AT, A, A, ux, — «2, -V3] 

These discrete symmetries imply that without loss of generality we can restrict the variables A > 0, and 
A > 0, since the dynamics in the other regions can be obtained by simply applying one or more of the 
maps above. The third and fourth symmetries listed imply that one can add additional constraints on the 
variables £x2,£i3,i'2 or v 3 ; however, in general there is no natural way to restrict any one of the variables, 
and hence we will not do so here. Note that any equilibrium point in the region v 2 > has a matching 
equilibrium point in the region v 2 < 0. 

2.2. Invariant sets. In this analysis we will be concerned with the following invariant sets: 

(1) T(VI h ): The general tilted type VL h model: |A| < 1. 

(2) Tx(FIh): A one-tilted type VL h model: |A| < 1, v 2 = v 3 = £ 12 = £ 13 = 0. 
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Table 2. The dimensions of the invariant sets for h ^ —1/9. The group parameter h is 
regarded as fixed. 



Dim 


2 


3 


4 


5 


6 


7 




B (VI h ) 


Ti, (VI h ) 


B{VI h ) 
T£{VI h ) 






T(VI h ) 



(3) Tifi(VIh)- A one-tilted diagonal type VIh model: v 2 = v 3 = £ i2 = £13 = £23 = A = 0. 

(4) N^iVh): A class of tilted type VI h models with W° = (sec cq. pT23|) for definition): |A| < 1, 
N ab v a v b = (/i^ -1/9). 

(5) T^(VIh)'- A two-tilted type VI^ model. This is the fixed-point-set of <^j~ and is given by V3 = £13 = 
E 23 = A = 0. 

(6) T 2 ~{VIh)'- A two-tilted type VIh model. This is the fixed-point-set of <p^ and is given by v 2 — £12 = 
£23 = A = 0. 

(7) B(VI h ): Non-tilted type VI A : |A| < 1, v ± = v 2 = v 3 = £12 = £13 - 0. 

(8) B a (VI h ): A class of diagonal non-tilted type VI h models (n a a = 0): V = £ i2 = £13 = £ 23 = A = 

(9) T{II): The general type II model: A = ±1, A = 0. 

(10) B{I); Type I: N = A = V = 0. 

(11) 5T(J): "Tilted" vacuum type I: O = JV = A = 0. 

Regarding N ± (VIh) 1 this invariant set is not a manifold; it is similar to the light-cone in 2-dimensional 
Minkowski space. Therefore, 7V ± (X / // l ) — Ti(VIh) consists of 4 disconnected pieces. By the symmetry ^4, 
these are actually only two inequivalent pieces. Here, we choose N ± (VIh) such that 

T+{VI h ) c N+(VI h ), T^{VI h ) C N-(VI h ) 

Since N+(VI h ) n iV _ (^4) = T t {VI h ), both ^+(^4) and N~{VI h ) are invariant sets. 
We note that the closure of the set T(VIh) is given by 

(2.22) T{VI h ) = T(VI h )UT(II)U B(I)UdT{I). 

Since the boundaries may play an important role in the evolution of the dynamical system we must consider 
all of the sets in the decomposition (|2.22|) . 

2.3. The case h = —1/9. A few comments are in order when h = —1/9. Let us consider the constraint 
equations (|2.19[) and (|2 . 20[) as a linear map 

L: (£i 2 ,£i 3 ) h-> (v 2 ,v 3 )/G + , 

where L is considered given in terms of A, N, A and Q. For h 7^ —1/9, det(L) ^ and the image of L is 
2-dimensional. However, for h = —1/9, det(L) = and the image of L is 1-dimcnsional; hence, in this sense 
the constraint equations are degenerate. This implies that (^2,^3) has to be restricted to a 1-dimensional 
submanifold. We will therefore say that the general type VI_x/9 model only allows for 2 tilt degrees of 
freedom. 

On the same token, this also mean that (t>2, V3) — does not necessarily imply that (£12, £13) is zero. In 
particular, for the non-tilted models this implies that we may have an additional shear degree of freedom; 
these models have usually been called the exceptional case and are denoted with an asterisk: B(V T^yg). 

It is advantageous to redefine iV ± (y/; l ) for h = —1/9 because, as can be shown, N a i,v a v b is identically 
zero for h = —1/9. Let us instead define 

D = [A(£? 2 + £^ 3 ) + 2£ 12 £i 3 ] , 

and define N ± (VI^ 1/9) as the set of points where D — 0. This is the natural generalisation of N^iVIh) 
to h = —1/9. We also note that for the tilted models, there is an exceptional case of the one-tilted models 
Ti(VI_ 1/9 ) which could be denoted Ti(VT* 1/9 ). However, Ti(VF 1/g ) = N-(VI_ 1/g ) where iV~(V7_i/ 9 ) 
is defined above. Therefore, we keep the notation N~(VI_i/q). 
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Table 3. The dimensions of the invariant sets for h = —1/9 in terms of the number of 
tilt degrees of freedom. The left-most column indicates the specialization in terms of the 
G*2 cosmologies. Here, HO means hypersurface orthogonal and an asterisk indicates the 
exceptional case where the models aquire an addition shear degree of freedom. 



Dim 


2 


3 


4 


5 


6 


7 


Gen G 2 






vacuum* 


0* tilt 


1*, 2 tilt 


2 tilt 


HO KVF 


vacuum* 


0* tilt 


1* tilt 









2.4. Fluid Vorticity. The various invariant subspaces can also be categorised in terms of the (-Hfluid- 
normalised where -fffluid = (1/3)m a '. a1 ) fluid vorticity, W a . The vorticity of the fluid for the type Vlft, models 
is given by: 

(2.23) Wa = ^B ( Nabvb + £abcvbAC + y2 N ^ b v c vA , W = -v a W a , 



where 

B = 



l ~W 2 + V 2 S + 2A a v a ) 
Vl~V 2 [l-{j-l)V 2 } 
For the invariant sets: 

(1) T(VIh): General vortic type VI?,, where all components W a can be non-zero. 

(2) N±{VI h ): W° = W 1 = 0. 

(3) T+(VI h ): W° = W 1 = W 2 = 0. 

(4) T^{VI h ): W° — W 1 — W 3 = 0. 

(5) Ti(F7 h ): W° = W a = 0, non-vortic. 

(6) B(VI h ): W° = W a = 0, non-tilted and non-vortic. 

In the special case of h = —1/9 further components of the vorticity are zero due to the degeneracy of the 
constraint equations as explained above. 

3. Qualitative behaviour 

3.1. Monotone functions. There are a number of monotone functions in the state space of interest. For 
< 7 < 6/7, there exists a monotonically increasing function Z\ defined by 

(3.1) Z x = aft 1 " 7 , H , 

Z[ = [2(l- 7 ) g + (2-7)+ 7< S]Z 1 . 

To see that this is a monotonically increasing function, note first that |5| < 2E 28J. Then, using the 
constraint equation (|2. 1T|) and (1 — E) > (1 — S 2 )/2, we can write 



2(l- 7 ) g +(2- 7 )+ 7 5 



7 (l- 7 )(^ 2 + 3) ( 



(3.2) > (6 - 7 7 )S 2 + 2(1 - 7 )(|N X | 2 + A 2 ) + ^ ' — ft, 

which is strictly positive for 7 < 6/7. Thus for < 7 < 6/7, Z\ is monotonically increasing as claimed. 

Theorem 3.1. For < 7 < 6/7, all tilted Bianchi models (with f2 > 0, V < 1) of solvable type are 
asymptotically non-tilted at late times. 

Proof. Use of the monotonic function Z\. □ 

In fact, the result that these models are asymptotically non-tilted at late times is true for all 7 < 1; this 
follows from further analysis and numerics, but is not covered by this theorem. This theorem is, in fact, true 
for all ever-expanding Bianchi models, including ever-expanding class A models. 

An immediate result of the above is the following corollary: 
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Corollary 3.2 (Cosmic no-hair). For Q > 0, V < 1, and < 7 < 2/3 we have that 

lim = 1, lim V = 0. 

Proo/. See pQ. □ 
Let us define 

(3.3) L» = iV|v| 2 + Re(N* x v 2 ), 

for which 

D' = (q + 2T + 2A Vl )D. 

This implies that D = is an invariant subspace and defines N^iVIh)- Moreover, the following function is 
a monotone function in T(VIh), h 7^ —1/9: 

(3 - 4) Z " = (1-^-7)11' 

(3.5) Z 2 = (57-6)(3-2A«i)Z 2 . 

This function is monotonically decreasing for 7 < 6/5 and monotonically increasing for 6/5 < 7. 
For a monotone function which is also monotone for h = —1/9 we define 

~ N 3 G 2 

(3.6) D = (1 — A 2 ) [A(S 2 2 + £ 2 3 ) + 2£ 12 £ 13 ] 

Using the constraint equations we have D = — (1 + 9/i)D so using D instead of D in Z2 gives the same 
equation of motion. 

In the subspaces Ti(VIh) and iV _ (y/_i/g) we have the monotone function: 

d 7) 7 = ^ Q 

1 J 3 A2G+(l-7 2 )i( 2 ^)' 

(3.8) Z 3 = -(2-7)(3-2i4«i)Z 3 . 

This function is monotonically decreasing in Ti(VIh)- 
The monotonic function Z% immediately implies: 

Theorem 3.3 (Future behaviour in T x {VI h ) and iV"(F7_i/ 9 )). For 2/3 < 7 < 2, Q > 0, A > 0, v\ < 1, 

V2 = V3 — we have that: 

either lim O = 0, or lim V = 0. 

r — >oo r — >oo 

This implies that all non-vortic type VI/t universes are either asymptotically vacuum or non-tilted at late 
times. 

In the remains of the section we will present all the equilibrium points including their eigenvalues, and 
discuss their stability. 

3.2. Equilibrium points. The Bianchi type VI/j models possess a wealth of equilibrium points. A full 
catalogue of the equilibrium points has up until now been lacking. Therefore, several of the equilibrium 
points we list are new. Due to the function Z2, the type VI ^ equilibrium points can be divided into 4 cases, 
according to whether D = 0, = (vacuum), 7 = 6/5, or V = 1 (extremely tilted). The equilibrium points 
in dT(I) are all unstable and will not be given here. 

3.2.1. B(I): Equilibrium points of Bianchi type I. 

(1) 1(1): S + = E_ = £12 = £13 = £23 = A= N = V = and = 1. Here, |A| < 1 and is an 
unphysical parameter. This represents the flat Friedman-Lemaitre model. 
Eigenvalues: 



2 The list in Apostolopoulos' paper |34| is incomplete. 
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3.2.2. T(II): Equilibrium points of Bianchi type II. All of the tilted equilibrium points come in pairs. These 
represent identical solutions (they differ by a frame rotation); however, since their embeddings in the full 
state space are inequivalent, two of their eigenvalues are different. All of these equilibrium points have an 
unstable direction with eigenvalue — 2\/3£23 corresponding to the variable A. These equilibrium points are 
given in [22] . 



3.2.3. T(VI h ): D = case (ft > 0, V < I). 



(1) C(h): Collins perfect fluid solutions, h< 0, 2/3 < 7 < ={ 
£ 12 = E 13 - S 23 = A = V = 0, £+ = -i(3 7 - 2), £_ 
A = V^ShN, = f [(2 - 7) + /i(3 7 - 2)]. 

This equilibrium point is in B(VIu). 

Eigenvalues: 



2(1-/1) 



-3h 



-3h 



(37-2), N- 



2 — JL 

16 



(3 7 -2)(2- 7 ), 



Ai + A 2 = A 3 + A 4 = --(2 - 7), AiA 2 > 0, A 3 A 4 > 0, 

As = _(2-t)^ A 6 , 7 = ^[(5 7 -6)±v / ^(37-2)] 
7 2 7 



(2) K+(h) 

S13 = S 2 3 



2(3-\Z^7T) 
5-3 



?<7<|, -|<^<0 



and 



2(3-v^ft) 2(1+3^7?) 



k = V^h: 



A 0. £ . - _ . 



where 



A = 576 k 2 7 2 (7 - 1) (1 + 3 fc) 2 (-36 k 2 - 24 fc - 4 + 120 fc 2 7 + 160 fc 7 + 40 7 - 133 fc 2 7 2 

- 182 fc 7 2 - 37 7 2 + 49 k 2 7 3 + 54 k 7 3 + 9 7 3 )(6 - 5 7 + 3 fc 7 - 2 /fc) 2 



a = 16 (35 7 - 36) (7 - 2) + (-2880 7 3 + 12192 7 2 - 15936 7 + 6528) k+ 
(33792 7 2 + 8064 - 18432 7 3 - 26880 7 + 3600 7 4 ) k 2 



288 (7-I) (15 7 3 - 47 7 2 + 46 7 - 12) k 3 + 432 7 (3 7 - 2) (7 - l) 2 A: 4 



23616 7 2 ) k 



b = 8 (7 - 2) (35 7 - 36) (3 7 - 2) + (-6528 + 11976 7 3 + 20256 7 - 2160 7 4 
+ (14280 7 3 + 25248 7 - 8064 - 2544 7 4 - 28800 7 2 ) k 2 + 
(-10248 7 3 + 19392 7 2 + 1584 7 4 - 14112 7 + 3456) k 3 + 144 7 2 (3 7 - 2) (7 - 1) k 4 

V _ (3 7 -2)+4S + + fc(l-2S + )(5 7 -6) y2 (l-2S + )(C +CiS + ) wWp 

Co = (2 -7)(57-4)(3 7 - 2) + 2(-197 3 + 607 2 - 647 + 24)A:-7(37-2)(57-6)A: 2 



Ci 



4(2 - 7) (57 - 4) + 8(5 7 3 - 24 7 2 + 32 7 - 12)fc - 4 7 (3 7 - 2)(2 7 - 3)k 2 



N 2 - [(37-2)+4S + ][4(57-6)(3fc(7-i)-i)s+-(37-2)(5 7 -6)+3fc(77-6)(2-7)] a - JZu N T he exnressions for 
i\ — i2fc a 7 a [(37-2)/c-(5 7 -6)] ' 71 ~~ VoKiv. i ne expressions ior 

ft and V 2 are given in the appendix, eqs. (|A.4|) and (|A.5[) . 

This equilibrium point lies in T% (Vlh)- 

Eigenvalues: 



Re(Ai, 2 ,3,4) < 0, A 5 +A 6 = -2(l + £+), A 5 A 6 > 0, 



At = -(5 7 -6)(3-2A Wl ). 



This equilibrium point is stable in N + (VIh), but always has one unstable eigenvalue in T(VIh). 



(3) K-{h): 



5+3%/^Ti 



< 7 < 



4 < h< 



and 



1 < -7 < — i— 



h) 



-Kh< 



-h: 



This point is given via the expressions for TZ + (h) above with k < but £+ = — an d using the 
symmetries fa, 4>i and 04 (thus interchanging Eia and E13, and ^2 and W3 so that Ei 2 = v 2 = 0). 
This equilibrium point lies in T^~(VIh)- 
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Eigenvalues: 

Re(Ai, 2 , 3 ) < 0, Rc(A 4 ) = 



< 0, -l/9<ft<0 
>0, -K/K-l/9, 



A 5 + A 6 = -2(1 + E+), A 5 A 6 >0, A 7 = i(5 7 -6)(3-2A Wl ). 

This equilibrium point is stable in N~(VIh) for -1/9 < ft < 0, is stable in T(VI h ) for 2 ^+^ } < 
7 < |, — | < ft < and has one unstable eigenvalue otherwise. The eigenvalue A4 corresponds to a 
direction in T~(VIh), while A7 corresponds to the quantity D, defined earlier. 

3.2.4. T{VIh): Vacuum case (Q = 0). All of these equilibrium points are plane wave solutions and exist for 
< 7 < 2, ft < 0. Moreover, they all have D = 0, and 



n = £12 - E13 = £23 = 0, £_ = N = v/-£+(l + £ + ), A = (1 + £+), -!<£+< 0, |A| < 1. 



The group parameter is given by 3ft£+(l — A 2 ) = (1 + £+). It is also advantageous to introduce r = \f\ — A 2 
and the parameter K 2 = — 1/ft (K can have either sign). This implies that we can write 



We will also define p by 



p = v\+v\. 



For all equilibrium points with = 0, three of the eigenvalues are always 

Ai=0, A 2 ,3 = -2[(l + S+)±2iV3AiV]. 
The equilibrium points are then determined by the tilt velocities: 

(1) C(h): vi = V2 = V3 = 0. These represent 'non-tilted' plane waves and lie in B(VIh)- 

Eigenvalues: 

A 4 = -k^T, b(* 2 + **) - 2(^ + r 2 )] , A 5 = [2{K 2 + 2r 2 ) _ , {K i + 3r2) ] , 

K7 = - K 2 3 +3r 2 + ^ ± \ K \ r2 ) 7(^ 2 + 3r 2 )] 

(2) C(h): vi = 7(g2+3 2 r r ' ) ( ; 2 _ ( i y +2r2) , v 2 = v 3 = 0, ^g'+ffi < 7 < 2. These represent 'intermediately 
tilted' plane waves and lie in T\(VIh)- 

Eigenvalues: 

3(if 2 + r 2 )(2- 7 ) 



A 4 
A 5 = - 

^6,7 = 



(iP + 3r 2 )( 7 - 1)' 

3(if 2 + r 2 )[{K 2 + 5r 2 ) 7 - 2{K 2 + 3r 2 )][(if 2 + 3r 2 ) 7 - 2(K 2 + 2r 2 )] 
(K 2 + 3r 2 )( 7 - \)[{K 2 + 3r 2 ) 2 7 - 2(/f 4 + 4K 2 r 2 + 5r 4 )] 
3(X 2 + r 2 )[37-4± \K\(2-j)} 
2(^ 2 + 3r 2 )( 7 - 1) 



(3) C±(h): v\ = ±1, v 2 = vs = 0. These represent 'extremely tilted' plane waves and lie in T\(VIh)- 
Eigenvalues: 

6(K 2 +r 2 ) 

C + (h): A 4 = 0, A5 = 2A 6 = 2A 7 = 7 i^-J 



C-{h) : A 4 — - - A 5 



(K 2 + 3r 2 

12r 2 6[(X 2 + 5r 2 ) 7 -2(^ 2 + 3r 2 ) 



if 2 + 3r 2 ' & (if 2 + 3r 2 )(2- 7 ) 

^--^3^ (r2±2|X|r2 "^ 2) - 



10 



S HERVIK, R J VAN DEN HOOGEN, W C LIM AND A A COLEY 



(4) T+{h): Here K > and 

(K 2 + r 2 ) [2(2 + K)-(3 + K)j] [ 7 (if 2 + 3r 2 ) - (K 2 + Kr 2 + 4r 2 ) 



1 (K 2 + 3r 2 )-{K 2 + Kr 2 + Ar 2 ) 2 2 
vi = r2(3 _ 2l + K) . v 2 -v 3 =si g n(K)pr 



P 

where 



r 4 (l + K){3-2~f + K) 2 
2(K + 2) ^K 2 + {A + K)r 2 „ 2 

£+^< 7 <™ for r 2 <_K<r { r + ^ 

^ + {A + K)r 2 K 2 + 3r 2 ^— 

K 2 + 3r 2 - 1 ~ K 2 + (2-K)r 2 ' V+V 

These represent 'intermediately tilted' plane waves and lie in N + {VIh) (for A = they lie in 
T 2 + (VI h )). 
Eigenvalues: 

3(K 2 +r 2 )F(K,r, 7 ) 

A6 + A7 



A6A7 — 



(3 - 2 7 + K) (K 2 + 3r 2 )G(K, r, 7) ' 
18(X 2 + r 2 ) 2 [2(2 + K) - (3 + Jf) 7 ] 



(X 2 + 3r 2 ) 2 (3 - 2 7 + K)G{K, r, 7) 
x [7(,ft: 2 + 3r 2 ) - (X 2 + XV 2 + 4r 2 )] [(K 2 + 3r 2 ) - j(K 2 + (2 - JQr 2 )] , 

where we have defined 

G(K, r, 7) = (-K 4 + 3r 4 - AK 2 r 2 + 2K 2 r 4 + 8r 4 K) 

+ (7K 2 r 2 - 9r 4 K - K 2 r 4 + 2K 4 - K 3 r 2 - 2r 4 ) 7 + K(-K + r 2 )(K 2 + 3r 2 ) 7 2 
F(K, r,-y) = 2(K + 3) (-K 4 - 5K 2 r 2 + 2K 2 r 4 + 6r 4 K - 2r 4 ) 

+ (16K 3 r 2 - 65r 4 K + 17K 4 + 78K 2 r 2 - 38K 2 r 4 + 5K 5 + Air 4 - Ar 4 K 3 - 2ifV) 7 

+ (-40r 4 - &JK 2 r 2 - 16K 4 + r 4 K 3 + h\r 4 K + 3K 4 r 2 + 20K 2 r 4 - 4K 5 - 4if 3 r 2 ) 7 2 

- (K 2 + 3r 2 ){K 2 r 2 + 5Kr 2 - 4r 2 - 5K 2 - K 3 )-/ 3 

(5) !F~{h): These are given by the same expressions as for F + (h) but with K < 0. The range of 7 is as 
follows: 

K 2 + (4 + K)r 2 K 2 + 3r 2 



K 2 + 3r 2 K 2 + (2-K)r 2 > ^ ^ < ^ - vWl), 

^^^<^, - K r-^)<^<0, 

These represent 'intermediately tilted' plane waves and lie in N~(VIh) (for A = they lie in 
T 2 -(VI h )). 

Eigenvalues: As for T + {h) but with K < 0. 
(6) £+{h): Here K > r(r + Vr^+l) and 

r 2 (l + K) 2 2 
Vl = ~ K{K-r 2 Y v 2- v 3=^MK) P r 

2 _ (K 2 +r 2 )(K 2 -2r 2 K -r 2 ) 
P v i~ K 2 (K-r 2 ) 2 

These represent 'extremely tilted' plane waves and lie in N + {VIh) (for A = they lie in T^iVIh))- 
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Eigenvalues: 



A 4 = K(K-r 2 )(K 2 + 3r 2 ) ' 5 = 2%/3 ( 1 -vi)sign(K)Nr, 

3(K 2 + r 2 )[{5r 2 +K 2 ) 1 -4r 2 (2 + K)} 

A(5 + A7 — — - 



A fi A 



6^7 



K(K-r 2 )(K 2 + 3r 2 )(2- 7 ) 
18(K 2 + r 2 ) 2 (if 2 - 2if> 2 - r 2 )[(if 2 + 2r 2 - Xr 2 ) 7 - (if 2 + 3r 2 
if 2 (if - r 2 ) 2 (Jf 2 + 3r 2 ) 2 (2 - 7) 



(7) £p{h): These are given by the same expressions as for £p(h) but with K < r(r — \/r 2 + 1). 

These represent 'extremely tilted' plane waves and lie in N~(VIh) (for A = they lie in T^(VIh))- 
Eigenvalues: As for £p(h) but with K < r(r — y/r 2 + 1). 

2.5. T(VI h ): 7 = 6/5 case (ft > 0, V < 1, D ^ 0). 

(1) B(h), -1/9 < h < 0, 7 = 6/ 5, fc = V ^: 

E 23 = 0, A = x/ZkrN, r = y/l - A 2 , u x = \A, £_ = |iVA, £+ = -§(1 - A 2 ), 

n = T [3(1 - A 2 ) — TV 2 ] , (Eia.Eia) = E x (contain VO, (va,T> 3 ) = v ± (cos<f>,sin(p), 

E» = -\V 2 {1 - A 2 ) + i(2 + y 2 )iV 2 - ^(1 - A 2 ) 2 - lN 2 A 2 , vl = V 2 - §A 2 , and $ is given by 

-7r/2 < ip < and 

, , N , -12MiV 2 + l6MN 4 k 2 r 2 - 15M 2 v 2 , + 15E^7V 2 r 2 (9fc 2 - 1) 
cos(2-0) = kr 



55^ [2A/ + iV 2 r 2 (l-9fc 2 )] 



where M = 1(1 - A 2 ) - 2iV 2 . 



The tilt velocity V and the angle <fi are given in terms of (fc, N, A) (see Appendix, eqs. (|A.1[) and 
(|A.2[) ^) . The variable A is determined implicitly by a sixth-order polynomial equation P(k, N, A) = 
in A 2 where N and k appear as parameters. This polynomial is also given in the appendix, eq. (|A.3p . 
There is a unique solution in terms of A, up to symmetries, as long as the range of the parameters 
fc and N are as follows: 

< fc < I, N 2 ^. <N 2 < 12 



3' -lower--' -25-9/c 2 
where 

2 _ 65 - 60fc + 27fc 2 - 5(1 - 3fc)V9fc 2 + 6fc + 49 
lowcr ~ (5 + 3fc)(-18fc 3 - 15fc 2 - 20fc + 25) ' 

If and only if the parameters fc and N obey the bounds above, there exists an equilibrium point 
for the system. Hence, for every — 1/9 < h < 0, 7 = 6/5, there exists a one-parameter family of 
equilibrium points, parameterised by N. In FigQ]we have plotted the surface P(k, N, A) = defining 
the solutions B{h). 

This family of non- isolated equilibrium points connects lZ~(h) for —1/9 < h < and 82(h) at 
7 = 6/5. 

Eigenvalues: Numerical calculation of the eigenvalues shows: 

Ai = 0, Re(A 2 , 3 ,4 : 5,6,7) < 0. 
Hence, this line of equilibria is stable in T(VIh) whenever it exists. 

2.6. T(VI h ): Extremely tilted case (Q, > 0). 

(1) £i(h), -1/9 < /i < and < 7 < 2, k = \/^h: 

v _ v _ _ n V - (15fc 2 -2fc-5) v _ v/3(21fc 2 -10fc-3) v _ V3-7fc 2 (l-3fc) 
^13 - ^23 - V 3 - U, lu+ - - 4 (l+fe)(3fe-2)' ^~ ~ 12(l+fc )(3fc-2) ' ^ ~ yi2(2 -3fc)(l+fc) ' 

N = \l ^3^XT W, A = V3kN, Vl = —\ .J^fi iv, v 2 = J rTo^tS t and n _ " l 



2(2-3fc)(fc+l) 2 ' ^ — V Oh-iV , t/i — y 6(2-3fe)(l-fe) > V2 ~ \ 6(2-3fe)(l-fc) <mu " — 2(2-3fe)(l + fe) 2 ' 

This equilibrium point lies in T^(VI h ). 
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Figure 1. The equilibrium points B(h): Here the surface P(k,N,X) = is plotted in 
(k, N, r)-space (recall that r = yl — A 2 ). 




Eigenvalues: 

(3 - 7k 2 ) 



Ai = - 



A 6 = - 



2(2-fc-3fc 2 ) 

(3-7fc 2 )(2 7 -3) 
(2-3fc)(l + fc)(2- 7 ) : 



^2,3 — ^45 — 



4,5 



3(1 — k) 2 



A 7 = 



4(2-3fc)(l + fc) 

3(3 - 7k 2 ) 
2(2-fc-3/c 2 ) 



1±4 1- 



32(3-7fc 2 )(l-3fc)(2-3fc) 
9(1 -k) 4 



This equilibrium point is stable in N + (VIh) for 7 > 3/2, and always unstable in T(VIh)- 



(2) £f(ft), -1/9 < ft < and < 7 < 2, k 



-ft: 



This point is given via the expressions for £^(h) above with k < and using the symmetries <j>2, <Pi 
and 4>4 (thus interchanging E 12 and E 13 , and v 2 and v 3 so that S 12 = v 2 = 0). 

This equilibrium point lies in T 2 (VIh)- 
Eigenvalues: Expressions same as for £^~(h) but with k < 0. 
Stable in N~(VIh) for 7 > 3/2, always unstable in T(VI h ). 
(3) £ 2 (h), -1/9 < ft < and < 7 < 2: 

y/-3hN, «i 



E 23 = A = 0, E+ = - 



iV 



_ 6(l+9/Q 
52 ~~ 25+9/i 



S_ = |v^3ftiV 2 , A 



25+9h ' Y 25+9/i' 3 V ' oi — 3 

, (£12, £13) = Ej_(cos^,sin^), (v 2 ,v 3 ) = vj_(cos <p, sin </>) where = 15(1 (25 ^ ( 9 5 v , 2 27fe) , 

v l = T5+W' cos2 ^ = 27 ^27h h) ' */2<iP< tt, and cos20 = < < tt/2. 

Eigenvalues: 



2 v /T 3ftX 

0(5-s 
(25+9/1) 2 



Al,2 



A 



5,6 



^3,4 — 
15(1 



3(5 -3ft) 
" (25 + 9ft) 

•ft) 



1± 



20(1 + ft)(25 + 9ft) 



(25 + 9ft) 

This equilibrium point is stable for 7 > 6/5. 
(4) £ 3 (ft), -1 < ft < -1/9 and < 7 < 2, k = -\ 

S12 - S23 = A = 0, £+ = - % + fc 6fc 1 ^ 1) , E 

v /3(l-4fc-fc 2 )(l+fc) 2 _ 3(l-/c 2 )(l+3fc) 
^13 - V 4(l-3/c) 2 ' " - 2(l-3fc) 2 

This equilibrium point lies in T 2 (VIh)- 



A 7 



(5 -3ft) 2 

15(1 + ft) 
25 + 9ft 



(5 7 - 6) 



V3 (l+6fc+fc 2 ) 
4 1-3/s • 



3(l-fc) 



2(l-3fc) 



7: 



i4 = V^hN, 



_ (i+fc) 



, ui = - 



w 2 



— / 1—4/ 

- v~ 



4fc-fc 2 



— 
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Eigenvalues: 
Ai - 

A 4 = 



3(l-4fc-fc 2 ) 

2(1 - 3k) ' 2 ' 3 

3(l-4fc-/c 2 )( 7 -l) 



(l-3fc)(2- 7 ) 

This point is stable for 7 > 1. 
(5) TlC{h), < 7 < 2, h< 0: 



3(1 -k) 2 
4(1 - 3k) 

A5,6 = A2.3, 



1± 



A7 = Ai 



32(l-4fc-/c 2 )(l + fc) 



£ 12 = £13 = £23 = v 2 = v 3 = 0, £_ = N = (1 - ^) A /-E + (1 + E+), A = (1 + £+), «i = 1, 
-1< E + < 0, < £ < 1, |A| < 1. The group parameter is given by 3/iE+(l - A 2 )(l - £) = (1 + £+). 
These equilibrium points correspond to a non-vacuum plane wave with an extremely tilted fluid. 
These equilibrium points lie in Ti(VIh) and are always unstable. 

3.3. Special equilibrium points for h = —1/9. 

(1) CT: The Collinson- French (Robinson- Trautmann) solution is given by: 



£+ : 
£12 



-^23 = = 



1 



" 3>/3' 
A = 0. 



^13 



1 

71' 



.4 



1 

7!' 



As for the plane waves the w^-equations decouple and we can again treat these separately. The 
special case, h = —1/9 leads to the exact vanishing of one of the constraint equations, and hence the 
tilt-velocity can only have 2 independent components (and we set V3 = 0). 

There are the following equilibrium points associated with the Collinson-French solution: 

(a) CT : Vi = v 2 = 0, < 7 < 2. 

(b) CT 1± : Vl = -^^, 

^(9 7 -14) 
6(7^1) ' V2 



± V5(3 7 -4) (3-2 2 ) ± 3 

U2 ^ n/2(3- 7 ) ' 3 ^ 7 ^ 2- 

n 24- V6 < < 24+V6 
u ' 15 ^ ' v 15 ' 



(c) CT 2 : Vi 

(d) £C]F ± : ui = ±1, v 2 = 0, < 7 < 2. 
These equilibrium points, and their stability, were studied in [lj. 

(2) W: Wainwright 7 = 10/9 solution: 



£+ = 



1 



: A 



1 



3V3' 

y = 0. 



< £13 < 



15 



TV = 



54£ 2 3 , 



.4 



1 

71 



iV. 



3' 

£12 = £23 

This line-bifurcation was studied in [16] and found to be stable in T(VIh) for h = — 

3.4. Special equilibrium points for h = — 1. 

(1) ViJII); There is a special line-bifurcation of 'tilted' vacuum solutions for h = —1, 7 

■K E_=iV=^, A = -, 
4' 4 ' 4' 



1/9. 



£+ = — , E_ = N 
£12 = £13 = £23 = A = 0, vi =v 2 = 0, < u 3 < 1 
The extreme limit of this equilibrium point, lim^ 3 _,i Villi) = £pih) 



We will call 



£-(h) 



(h,r)=(-l,l) 



(h,r)=(-l,l) 



for £ illl) for simplicity. We will also define Vo(IH) = linxyg^o Villi) = 



Lih)\( h r ) = (_i i)- These solutions are, in fact, locally rotationally symmetric (LRS). 

Regarding the eigenvalues, there are 3 eigenvalues which are zero and the rest all have a negative 
real part. In order to resolve the stability properties of these equilibrium points one has to resort to 
centre manifold theory. 

(2) Voilll) = £(>OI(fc,r)=(-UV 

(3) £ -iIII)= £~ ih) 

(h,r)=(-l,l) 
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Figure 2. The regions of stability of the plane- wave solutions in the invariant subspaces 
T*2 (Vlh)- Here K 2 = — 1/h and K < corresponds to T^~(VIh) while K > corresponds to 
T^iyih). The dashed line corresponds to the lower boundary of the region marked C-(h). 
The loophole is shown for K > but has been suppressed for K < 0. 




-10 12 3 



3.5. New self-similar solutions. As we have provided a complete catalogue of equilibrium points for the 
tilted type VI/, model, it is of interest to discuss which of these solutions correspond to new solutions. The 
physically interesting solutions are the intermediately tilted non-vacuum points which correspond to exact 
vortic self-similar solutions of type VI?,. Self-similar solutions of type VI?, were studied by Rosquist and 
Jantzen [35] : however, only some type VLj solutions were found explicitly. The solutions 7Z~(h) appear to 
have been found by Apostolopoulos in [36]. Moreover, the h = —1/9 case of TZ + (h) appears also to be found 
by Apostolopoulos in [34]. However, the remaining solutions 7Z + (h) for h ^ —1/9 seem to be new and have 
not been given previously in the literature. Furthermore, the exact 7 = 6/5 solutions corresponding to B(h) 
also appear to be new (although the case h — was given in [3"Tl |2"8] ) . 

Regarding the extremely tilted non- vacuum equilibrium points, £ x (h), 82(h) and 83(h), these also corre- 
spond to exact solutions. We see these equilibrium points exist for all 7 € (0,2); however for 7 ^ 4/3 the 
physical interpretation of these is unclear. For 7 — 4/3 we can interpret these solutions as Bianchi type VI?, 
cosmologies containing null radiation. 

The remaining solutions listed are either known or implicitly known [21 [1] . 

4. Late-time behaviour 

Let us discuss the different late-time behaviours. 

4.1. The Realm of the plane-waves. For the type VII?, and the type IV universes, the vacuum plane- 
waves played a vital role in the future evolution [32] . These plane- wave solutions were also shown to possess 
a wealth of interesting phenomena, like attracting closed orbits and attracting tori. Since the type IV model 
is the h — > —00 limit of type VI?, models, we expect to find that the plane wave solutions play an important 
role for, at least, some of the type VI?, models. 

In the Bianchi type VI?, case all of the plane-wave equilibrium points reside in the invariant subspaces 
iV^VT?,). It is illustrative to consider the subspaces T 2 (Vlh) C N ± (VIh) (partly because it is easier to 
interpret; e.g., the figures for N ± (VIh) would be 3-dimensional) . The equilibrium points in T^(VIh) are 
those with r = 1. 
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Figure 3. Magnified region of Fig[2] showing the loophole in T 2 (VIh)- The curve marked 
F = is defined by F(K,r,j) — 0, where the function F is defined in section 15. 2. 41 item 
(|4|). Part of this curve marks the threshold value of stability for the equilibrium point T~ (h). 




^ ^ x A" 

-0.42 -0.4 -0.38 



In Tf{VI h ) the regions where the various plane-wave equilibrium points are stable are depicted in Fig. 
[2J The dashed line indicates the lower bound of stability for We can see that there are only stable 

plane-waves for h < -1 in T^{VI h ) and h < -1/9 in T}{VI h ). For the fully tilted models (i.e., in T(VI h )), 
only the plane- waves in T^~(VIh) remain stable (i.e., the ones with — 1 < K < in Fig. [2]); all plane waves 
in T+(VI h ) are unstable in T(VI h ). 

4.1.1. Hopf-bifurcation and the loophole. The loophole is defined as the set of plane-waves where the variable 
7 satisfies the following: given K and r = Vl — A 2 where r(r — vV 2 + 1) < K < r(r + y/r 2 + 1), then 
7o < 7 < 2(K 2 + 3r 2 )/(K 2 + 5r 2 ) where 70 is implicitly defined by F(K,r,-fo) = 0, where F is defined in 
section |3. 2. 4[ item ((4]). In the loophole, there are no stable equilibrium points, so we need to look for other 
potential candidates for late-time attractors. 

In [U [22] we noted that the corresponding loopholes in type IV and type VII^ models had closed curves 
and tori as attractors; hence, we need to look for closed curves in the type VI/j loophole as well. The 
important observation regarding the nature of the attractor can be seen from the eigenvalues Xq and A7 for 
the equilibirum point T (h). As we vary 7, we note the following (sufficiently close to 70): 

• 7 < 7o: A 6 + A 7 < 0, A 6 A 7 > 0. 

• 7 = 70: A 6 + A 7 = 0, A 6 A 7 > 0. 

• 7 > 7o: A 6 + A 7 > 0, A 6 A 7 > 0. 

This implies that the real values of Ag, 7 change sign while the imaginary values remain non-zero; this is an 
indication of a possible Hopf-bifurcation. 

Consider a 2-dimensional dynamical system given by the complex variable Z. The normal form of a 
Hopf-bifurcation can be written 

(4.1) Z' = {\ + b\Z\ 2 )Z, 

where b is some complex number and A is a parameter. If Re(6) is negative then there are stable closed 
orbits for A > 0. In order to determine whether our system experiences a Hopf-bifurcation as we vary K, r 
and 7 we thus need to expand to 3rd order in the variables. This has proven to be difficult in practice due to 
the complicated dependence on the parameters K, r and 7; however, there is analytic as well as numerical 
evidence that the fixed points J^ ± (/i) do indeed experience a Hopf-bifurcation and produce a stable closed 
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orbit as 7 passes through 70 for the following ranges: 



T + (h): r 2 <K < r(r + y/r 2 + 1), 

F~{h): r{r - yjr 2 + 1) < K < 0. 

This does indicate that there are attracting closed curves even outside the loophole; however, outside the 
loophole these attracting curves will co-exist with attracting equilibrium points. Inside the loophole, on the 
other hand, only the closed curves can be attactors as all the equilibrium points are unstable. The family of 
such attracting closed orbits, regardless whether they are outside or inside the loophole, will be referred to 
as the Mussel attractor. For numerical plots of the Mussel attractor outside the loophole, see Figs. l8land[T4l 
Let us study these closed curves in more detail. As in [22| , we consider the tilt-equations in a plane- wave 
background. We utilize the identity 

(v 2 2 +v 2 3 ) 2 = (2v 2 v 3 ) 2 + (v 2 2 -v 2 3 ) 2 , 

and define (x,p,0) by 

(4.2) x ~ vi, p = v 2 +v 3 , 2v 2 Vs = psinO, (v 2 — v 3 ) = pcos9. 
Then we have for the reduced system 

x' = (T + 2S;)x - (A* + V3N* cos 0)p, 
p = 2(T A*x-acos9) p, 

(4.3) 9' = 2a(A* +sin0), 

where a = \/3(l — x)N* and, by use of the discrete symmetry 04, we can assume that < 9 < 2n. 
Furthermore, these variables are bounded by 

(4.4) 0<p, V 2 =x 2 + p<l, |A*|<1. 

An asterisk has been added to the variables to emphasize that these should be thought of as the limit values 
for the full system. 

For a periodic orbit, c(r), with period T n , we introduce the average of a variable B: 

(4.5) {B) = ^-j>Bdr. 

We can also say something about the stability of a closed periodic orbit. For example, consider the evolution 
equation for f2, which we write as fi' = Afif2. Assume that c(r) is a closed periodic orbit with period T„. 
Then (An) indicates the stability of the closed curve with respect to the variable tt. 
We note that N^iVIh) intersect at sin6> = —A*, where 

N+(VI h ): cos6» = + v /l-(A*) 2 = +r, 



• N-(VI h ): cos6» = - (A*) 2 = -v. 

It is in these subspaces that the closed curves must exist. We also note that the variable 9 induces transitions 
N+(VI h ) -> N-(VI h ) which means that in T(VI h ), N+(VI h ) is unstable. 

Theorem 4.1. Assume that there is a closed properly periodic orbit, Cm parameterised by c(r), for the 
dynamical system Then 

_ 1 {K 2 + 3r 2 ) - (K 2 + Kr 2 + 4r 2 ) _ 3(K 2 + r 2 )[(5 + K)y - 2(3 + K)} 

{ ' {X) ~ r 2 (3~2 1 + K) ' (iP + 3r 2 )(3-2 7 + X) 

Proof. The proof is analogous to that in the type IV case; see [22] • □ 

Inside the loophole existence can also be proven: 

Theorem 4.2 (Mussel attractor). For A* < 1 and every K and 7 taking values in the type Vlh loophole 
there exists a closed periodic orbit, Cm, for the dynamical system 



Proof. The proof is analogous to that in the type IV case [55] and goes as follows. Inside the loophole there 
are no attracting equilibrium points. By restricting to the 2-dimcnsional set cos^ = ±r, we get a picture 
similar to FigHJ The equilibrium points (labelled A, B, C and F) are all hyperbolic and are sources/saddles 
as indicated. We note that the boundary consists of heteroclinic orbits (also drawn). We now see that there 
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Figure 4. Proof of the existence of the Mussel attractor. 



A 




B 



C 



is an orbit from B which cannot end at an equilibrium point and consequently must approach a closed curve. 
More explicitly, a future trapping region can be constructed as follows. We can use the shadowing theorem 
(see Proposition 4.2 in [5]) and the fact that the equilibrium points are all hyperbolic to show that there is 
an orbit, ci, originating from A which can be chosen to come arbitrary close to C and B (also shown in Fig. 
|4|). Similarly, there is an orbit, C2, originating from A and coming arbitrary close to B. We can then cut 
out neighbourhoods of the points A and B using the dotted curves, (1a, ds, so that ci, C2, (1a, cLb constitute 
the boundary of a future trapping region. Since the only equilibrium point inside this trapping region is 
the source F, the Poincarc-Bcndixson theorem implies the existence of a closed periodic orbit (in particular, 
there must be an orbit originating from x — p — that approaches a closed orbit). □ 

There is an important corollary that follows from this: 

Corollary 4.3. For h < — (3 — 2v2), there exist closed periodic orbits in T(VIh). 

The ramifications of this result is that considering the equilibrium points only is not sufficient to determine 
the asymptotic behaviour for these models. However, apart from the closed periodic orbits described above 
(the Mussel attractor), we have not found evidence for any other closed period orbits (for h ^ —1/9) 
important for the late-time behaviour. 

4.2. The Realm of non-vacuum vortic Universes. We note that all equilibrium points with £1 ^ in 
N ± (VIh) are, in fact, mT^(VI h ). The regions of stability of the equilibrium points in T^(VIh) are depicted 
in Fig. [5] For N ± (VIf l ), the picture is identical. 

In the fully tilted space T(VIh), we note that from the monotonic function Z2, if Q approaches a non-zero 
value at late times, then D — ► and V — > 1 for 7 < 6/5 and 7 > 6/5, respectively. Our analysis indicates 
that the future asymptote is non- vacuum for —1 < h < —1/9 and —1/9 < h < 0, for which we have the 
following behaviour (see FigJS]): 



(1) 



-Kh< -1/9: 
(a) § < 7 < 1: D -> 0, and 




V -> 0. Attractor: C(h). 

0, and V -> or V -> 1. Attractor: C(h) and 83(h). 
0, and V -> 1. Attractor £3 (ft). 



(2) 



->■ 0, and V -> 0. Attractor: C(h). 

— * 0, and V — > constant. Attractor: lZ~(h). 
and V — > constant. Attractor: 13(h). 
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(d) | < 7 < 2: D — > constant, and V — > 1. Attractor: £ 2 (ft)- 

4.3. The case h = — 1: Bianchi type III. The case /i = —1 is a special case which has to be treated 
separately. For 2/3 < 7 < 1 the Collins solution, C(h) is the only attractor. All eigenvalues have negative 
real parts and consequently the stability can be deduced from the linearised analysis. However, for 7 = 1 and 
7 > 1, the equilibrium points V(III) (including the non-tilted and extreme limits) and £~(III) have 3 and 
2 zero-eigenvalues, respectively. In order to determine the stability for these solutions one must therefore go 
to higher order. 

Consider first 7=1. Two of the zero-eigenvalues for V(III), < v 3 < 1 correspond to a non-trivial 
Jordan block: 




1 
' 



This means that the generic solution drifts along the line of equilibria. The solution drifts towards the V3 = 
solution, Vo (HI) ■ It is therefore of interest to study the centre manifold for V (III) . 

4.3.1. T2 (V I-i): The centre manifold. In our analysis of the centre manifold we restrict ourselves to the 
invariant subspace T 2 ~(VT_i). This subspace has only a two-dimensional centre manifold which makes the 
analysis more tractable and easier for illustrative purposes. It is useful to define the new variables 

£+ = i(-£+ + V3£_), 

(4.7) £_ = i(V3S+ + S_). 

The centre manifold can be found as follows. Define (x, y, z, w) by: 

(£+, N,v lt v 3 ) = (^+x,^-(l + y), z, w). 
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The variables f2 and £13 can be determined from the contraints. 

By introducing the variables (X, Y, Z, W) as follows: 

/ 3 3 3 3 

(X, Y,Z,W)= lx+-y,x- -y, -x - -y + w, -x + -y 

we can expand the differential equations for (X, Y, Z, W) to second order. The centre manifold can then be 
seen to be parameterised by (X, Z). Moreover, on the centre manifold (Y, W) = 0(3), where 0(3) means 
cubic in X and Z. The equations for (X, Z) to lowest order are then: 

X' = A 2 +0(3), 

(4.8) Z' = XZ + 0(3). 

Requiring f2 > 0, gives the approximate solutions at late times: 

t r 

where C is a constant. Substituting these back into the original variables, the decay-rates are found to be: 

1 



= 

2^3,0 



'13 



N = 
ft = 



V37-2' 



1 

r 

Wl = 



(4.9) w 3 



•^3,0 

r 

Here, the constant W3.0 is related to by C = 1 + W3.0. These decay rates have also been confirmed 
numerically. We have also simulated the system in the fully tilted space T{VIh) which does, indeed, confirm 
that the decay rates have the functional dependence as given above. The numerics also suggest that: 

(4.10) A = °(^)' ^ = 0(^3), v* = o(±y E»=o(^ 

Hence, it is plausible that: Vq(III) is the attractor for the fully tilted dust (7 = 1 ) type III model. 

For 1 < 7 < 2 the situation is easier. Again the the most promising candidate possesses zero eigenvalues; 
however, we note that 

£-{III)= lim £-(h)= lim £ 3 (h). 
p h->-l- p 

Also, for 1 < 7 < 2, £p(h), £z(h) are the attractors for h < — 1 and h > — 1 (sufficiently close to — 1), 
respectively; hence, we would anticipate that £~{III) is the attractor for h = — 1 and 1 < 7 < 2. Indeed, 
doing a centre manifold analysis for the subspace T 2 ~(W_i) we get: 

8+ 4 + (i), s-* + o(± 

(4.11) vl =o(^j, V3 =i+o^y n= °(l 

This confirms that the type III models with 1 < 7 < 2 at late times approach the extremely tilted vacuum 
solution £p (III) . 

Note that we have only presented here a centre manifold analysis for the invariant subspace T~(VI-\). A 
full analysis of these centre manifolds have been performed showing that these equilibrium points are indeed 
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the late-time attactors; however, due to the lengthy and complicated nature of the centre manifold analysis, 
this analysis will be presented elsewhere. 

There is an interesting connection with the type VIII models [35] worth pointing out. We note that the 
type III solutions approach an LRS universe at late times. In fact, type III LRS solutions also admit a type 
VIII action. Comparing the attractors for type III and type VIII we note there is a striking similarity: In 
terms of the metric and the matter (in a C° sense), the type III and type VIII attractors are the same. 
Explicitly, the curvature A 2 + N 2 in the type III model takes the role of —N\N in the type VIII model. In 
fact, by comparing the decay rates we even see that some of the decay rates are identical to lowest order. 
On the other hand, the two models do approach this attractor differently since the Weyl parameter diverges 
for type VIII while it is always bounded for type III. 

4.4. The case h = —1/9. Based on the eigenvalues of the equilibrium points, the equilibrium points are 
future attractors in the range specified pQ: 

(1) CT : 10/9 < 7 < 4/3. 

(2) £F 1± : §< 7 <§ + 4fI. 

(3) 0^2- Never an attractor. 

(4) ECT_: 24=?£ < 7 < 2. 

Including the analysis of the non-tilted equilibrium points we can thus conclude: 

• 2/3 < 7 < 10/9: Collins type VI_i/g perfect fluid solution is stable. 

• 7 = 10/9: Wainwright's type VP 1 ^ g , 7 = 10/9, solution is stable. 

• 10/9 < 7 < 2: The Collinson-French vacuum is stable. Moreover, for 4/3 < 7 the asymptotic tilt 
velocity is non-zero. 

It should be emphasised that the analysis of this model is not quite complete. Due to the exact vanishing 
of one of the contraint equations a different approach has to be taken regarding the numerical analysis. 
Furthermore, the question whether there are closed periodic orbits acting as attractors has not been addressed 
in this paper. Preliminary analysis indicates that this model also possesses closed periodic orbits 38 . 

5. Numerical Analysis 

5.1. Additional Numerical Support. In the qualitative analysis of the system of equations governing 
the evolution of the tilted Bianchi type VI^ models, it was found that there were numerous possible future 
asymptotic states. In addition, it was also found that in the invariant set T^iVIh) there is a portion of the 
parameter space, the loophole, in which there is no equilibrium point acting as a future asymptotic attractor. 
An extensive numerical analysis has been completed that further supports the qualitative analysis presented 
in the paper. In addition to obtaining numerical evidence for all behaviours for the general type VU, and in 
the special case of the invariant set T 2 f (Vlh), some very special and even somewhat surprising results have 
been found. Below we include a select group of numerical plots that illustrate some of the rich behaviour 
found in the type VI^ models. More details and numerical plots can be found in the companion article [39j . 

5.2. The invariant set T^iVIh)- Here, we are looking at the dynamics in the four-dimensional invariant 
set T+(VI h ). (Recall that E i3 = S 23 = A = v 3 = 0, and A = y/-3hN in the T^(VI h ) invariant set.) A 
six-dimensional set of differential equations for the variables {£+, £_, £12, N, v\, V2} was integrated using 
typical values for the parameters h and 7 that illustrate some different asymptotic behaviours. The constraint 
equations were used to check the accuracy of the numerical integrations. The built in differential equations 
solver ode45 in Matlab R2006a has been used to numerically solve all cases. The Relative Tolerance was set 
to 10~ 12 and the Absolute Tolerance was set to 10~ 20 or smaller. The constraint equations were used to 
determine the initial values of the tilt velocities i>i,«2. Two numerical plots, Figures [7] and [8] are included 
here that show some of the more surprising (albeit not totally unexpected) behaviour in the invariant set 
T^iVIu). Note the figures included here do not show typical behaviour for the entire range of parameter 
space for T^iVIh), but only illustrate what is typical for small neighborhoods around the selected parameter 
values. Indeed other behaviours are possible; there are additional situations when the future attractor is not 
unique, although, no others show evidence of closed orbits. For additional illustrations see the companion 
article [55] , 
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Figure 6. The asymptotic behaviour of tilted type VI?! universes. For h < — 1 we have used 
the solutions in T^iVIh) for illustrative purposes only (the figure should be 3-dimensional) . 
Again the loophole has been suppressed. 
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In Figure [7] we observe four periods of the variable V 2 and the corresponding Mussel attractor in the 
accompanying (vi,V2) phase portrait. The existence of this closed orbit was predicted, since the lack of 
a stable equilibrium point inside the loophole and the previous analysis of the Bianchi type VII^ and IV 
suggested that it must be true for at least some class of the type Vlh models. However, what was not predicted 
before the numerical analysis was completed was the existence of a closed orbit for parameter values outside 
the loophole. Figure [8] shows that it is indeed possible to have a Mussel attractor for parameter values outside 
the loophole. In addition, in this situation it can be seen that there are two attractors, the Mussel attractor 
which shows a periodic nature in the variable V 2 and the equilibrium point C-(h) which is extremely tilted. 

5.3. Fully tilted Bianchi VI^ space. Typical values of the parameters h and 7 are chosen to illustrate 
some of the interesting future asymptotic behaviour found in the fully tilted Bianchi VI/j models, see Figures 
I9lfl4l For the interested reader, additional figures can be found in the companion article [39]. In each 
numerical run we choose essentially the same set of initial conditions unless otherwise indicated. Given 
initial values for {£+, £_, £12, £13, £23, N, A}, initial values for A, V\,V2, «3 were determined by solving the 
constraint equations. For all cases except h = —5/9, the built in differential equations solver odeJf.5 in 
Matlab R2006a was used to numerically solve all cases. The Relative Tolerance was set to 10 -12 and the 
Absolute Tolerance was set to 10 -20 . The full 1 1-dimensional set of differential equations was integrated 
while the constraint equations were graphed (not included here) and used to determine when and if the 
numerical analysis breaks down. In the case h = —5/9, due to instabilities of the constraint surfaces, an 
alternative approach had to be taken. The built in differential equations solver odel5s in Matlab R2006a 
was used to numerically solve this case. Instead of integrating the 1 1-dimensional set of differential equations 
for the quantities {£+, £-, £12, £13, £23, N, A, A, Vx, V2, v%}, a 7-dimensional set of differential equations for 
{£+,£_, £12, £13, £23, N, A} and a set of algebraic constraints for A,vx,V2,V3 were integrated. Again the 
constraint equations were also graphed (not included here) and used to determine when and if the numerical 
analysis break down. The Relative Tolerance was set to 10~ 12 and the Absolute Tolerance was set to 10 -15 . 
In each phase portrait stars indicate the future attractor. 
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Figure[9]shows a partial picture of the non-isolated nature of the stable attractor B(h). In the (vt, \Jv\ + vf ) 
phase portrait, the lower bound (as seen in the phase portrait) of this non-isolated attractor is the equilib- 
rium point TZ~(h) (which is stable if the value of 7 is a little less than the critical value of 6/5) while the 
upper bound (not seen in the phase portrait) of this non-isolated attractor is the equilibrium point 82(h) 
(which is stable if the value of 7 is larger than the critical value of 6/5). Figure [TOl gives an illustration of 
a situation in which there are two stable equilibrium points in the fully tilted Bianchi Vlh models. Figure 
im shows the strongly attracting nature of the heteroclinic orbit between C(h) and 83(h). Figure fT2l not 
only shows the strongly attracting nature of the heteroclinic orbit between C(h) and F~(h), but also the 
non-isolated nature of the equilibrium point T~(h). We also note how the variables {£l,N, £ 2 ,A} obtain 
their equilibrium values much much faster than the tilt velocities. This behaviour is akin to the freezing in 
behaviour observed in the type Vllh and IV models [TJ H2] • Again we observe this same freezing in behaviour 
and the non-isolated nature of the equilibrium point S~(h) in Figure [TBI Figure [TTl illustrates something 
that was not predicted in the previous qualitative analysis; that is, the existence of two stable attractors, a 
stable equilibrium point £-(h) and a stable closed orbit. The parameter values h = —9,7 = 1.241 used in 
Figure [T4l would be considered outside the 3-dimensional analogue of the loophole. 

6. Discussion 

In this paper we have used dynamical systems methods and detailed numerical experimentation to inves- 
tigate the future asymptotic properties of SH Bianchi type Vlh cosmological models. We have determined 
all of the equilibrium points of the type VI^ state space. These equilibrium points correspond to exact self- 
similar solutions of the Einstein equations (some of which are new exact solutions) and play an important 
role in describing the general evolution of the system. The stability of all of these equilibrium points is 
also investigated. In particular, we have determined all possible late time behaviours for Bianchi type Vlh 
models. All of the possible future asymptotic behaviours for all Bianchi models (except Bianchi type IX 
models, which can recollapse) are now known and summarized in Table [H this Table is now complete. We 
note that the Bianchi type VI^ case is of special interest in that it is very complicated and contains many 
subcases, and many of the different future asymptotic behaviours found in previous work occur in these 
particular models. 

In particular, it was found that the vacuum plane-wave solutions play an important role in the future 
evolution of type Vlh models (as was the case in the type VII^ and the type IV universes [22] ). All of the 
plane-wave equilibrium points are found to occur in the invariant subspaces N^(VIh)- The regions where 
the various plane-wave equilibrium points were stable are depicted in Fig. [2j we recall that for the fully 
tilted models only the plane-waves in N^iyih) are stable. 

A loophole exists in which there are no stable equilibrium points, and hence it is necessary to seek 
other candidates for late-time attractors. It was noted that closed curves and tori act as attractors in the 
corresponding loopholes in type IV and type VII^ models [TJ [H] . Hence, we looked for closed curves in the 
type Vlh loophole. Calculations and numerical experimentation were found to provide plausible evidence 
that the system experiences a Hopf-bifurcation resulting in a stable closed orbit as 7 takes values in a 
range of parameter values (that includes values within the loophole, but also includes values just outside the 
loophole). This shows that there are attracting closed curves even outside the loophole; however, outside the 
loophole these attracting curves will co-exist with attracting equilibrium points and there are consequently 
a number of possible late time behaviours. Inside the loophole, however, only the closed curves can be 
attactors since all of the equilibrium points in the loophole are unstable. The family of attracting closed 
orbits (both outside and inside the loophole), are referred to as the Mussel attractor. 

In more detail, we analytically proved that in the type Vlh loophole there exists a closed periodic orbit (the 
Mussel attractor), Cm^), for the dynamical system. In addition, we provided convincing numerical evidence 
that closed orbits exist outside of loophole (see Figs. [7]and[H]). Clearly, an analysis of the equilibrium points 
alone is not sufficient for determining the future asymptotic behaviour in the Bianchi type VI^ models. 
However, we should also note that apart from the the Mussel attractor we have found no evidence for any 
other closed period orbits important for the late-time behaviour for h =/= —1/9. 

The Bianchi type III case (i.e., h = —1) is a special case which necessitates a separate treatment. For 
2/3 < 7 < 1 the Collins solution is the only attractor. However, for 7 > 1, the corresponding equilibrium 
points have multiple zero-eigenvalues, and a centre manifold analysis is required. We found that, in the 
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Bianchi 
Type 


Matter 


Future 
Attractor 


Comments 


T 

_L 


^ 7 z 


77 n 


Nn tilt allnwpH 


II 


|<7<f 

¥<7<f 
„, 14 
7= T 

if <7<2 


C«S(I7) 
W(77) 
L(II) 
£{II) 


Non-tilted Collins-Stewart 
Hewitt's tilted type II 
Tilted type II bifurcation 
Extremely tilted 


III 


!<7<1 

7-1 
1 < 7 < 2 


C(III) 
-Po(HI) 
8- (III) 


Non-tilted Collins type VI_i 

Non-tilted 

Extremely tilted 


IV 


1 <7<2 


Plane waves 


Tilted/non-tiltcd 


V 


1 <7<2 


Milne 


Tilted/non-tiltcd 


VL h 

h<-l 


2 ^ ^ 2(1-/1) 

3 ' l-3/i 

2(1 ~' l) < 7 < 2 


C(h) 
Plane waves 


Non-tilted Collins type VI ^ 
Tilted/non-tilted 


VL h 

-K/K-l/9 


!<7<1 

1< 7 < 2 ( 3 +^) 
2(3+v^) 2 


C(h) 

C(h) & £ 3 (h) 
£ 3 (h) 


Non-tilted Collins type VI,, 

Non-tilted or extremely tilted 
Extremely tilted, non-vacuum. 


vi_ 1/9 


1 < 7 < JU 

3 ^ ' 7" 9 

7=f 
^<7<2 


C(h) 
W 
CT 


Non-tilted Collins type Vl^/g 
Non-tilted Wainwright type VLjyg 
Collinson- French: Tilted/non-tilted 


-l/9</i<0 


2 2(3+v^) 

3 ^ ' v 

f^<7<! 
! < 7 < 2 


C(h) 

n-{h) 

B(h) 
£ 2 (h) 


Non-tilted Collins type VI^ 

Intermediately tilted 

Tilted type VI/j, 7 = 6/5 bifurcation 

Extremely tilted type VI^ 


VII ft 


l<7<2 


Plane waves 


Tilted/non-tiltcd 


VII 


2 < 7 < 4 

3 ^ ' ^ 3 

7=! 
I<7<2 


Pi{VIh) 
P 3 (VIIo) 
Pi{VIh) 


Non-tilted 

Radiation bifurcation 
Extremely tilted 


VIII 


|<7<1 
7=1 

1 < 7 < 2 


Pi{VIII) 

?2{VIII) 

Ei{VIII) 


Non-tilted 
Non-tilted 
Extremely tilted 



Table 4. The late-time behaviour of Bianchi models with a tilted 7-law perfect fluid (see 
the text for details and references). The case < 7 < 2/3 is covered by Corollary 13. 21 



invariant subspace T^iVIh) for 1 < 7 < 2, the type III models approach the extremely tilted vacuum 
solution £~(III) at late times. In addition, based upon a detailed centre manifold analysis in the invariant 
subspace T^iVIh) and numerical simulations in the fully tilted space T(VIh), we concluded that it is 
plausible that Vo(III) is the attractor for the fully tilted dust (7=1) type III model. 

As noted above, comprehensive numerical experiments were carried out to complement and confirm the 
analytical results. All of the numerics, which serve as confirmation of all of our analytical results, are 
summarized in a companion article [39) . In the present paper we have included several figures (from 39J), that 
provide typical results (confirming the analytical results) and some specific figures that contain interesting 
phenomena. 



We should also point out that the analysis of the special Bianchi type VI_]yg case is not complete either. Due to the exact 
vanishing of one of the contraint equations a different approach has to be taken regarding the numerical analysis. Furthermore, 
the question whether there are closed periodic orbits acting as attractors in this case has not been addressed in this paper; 
however, a prelimiary analysis indicates that such attracting closed curves do indeed exist |38| . 
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Appendix A. Nasty expressions 

A.l. Expressions determining V, </> and A for B{h). The following expressions are produced using 
Maple. 

(A.l) V 2 = 2 A 2 (-1728 - 724275 A 6 r 6 k 4 + 178848 A 4 k 4 r 4 - 11664 k 4 r 4 N 2 - 900 A 6 r 8 k 2 
+ 1350 A 8 r w k 4 + 26244 A 6 fc 8 r 8 + 32940 k 4 r 8 A 6 - 39690 k 6 r w A 8 
+ 13122 k 10 r w A 8 - 625 A 6 r 6 + 34560 A 2 - 31104 /c 2 r 2 N 2 + 150 7V 4 r 6 
+ 80000 N 6 + 368550 A 8 fc 6 r 8 - 525600 k 2 r 4 A 6 - 170640 A 4 k 2 r 2 
- 25000 k 2 r 6 N 8 + 414720 N 6 k 4 r 4 + 72075 N 6 k 2 r 6 - 77031 A 6 fc 6 r 6 - 100800 A 4 
+ 2592 k 2 r 2 + 864 r 2 - 90000 N 6 r 2 - 126846 A 8 k 8 r w - 8730 A 4 k 2 r 6 

+ 320400 k 4 r 6 A 8 - 27054 N 4 r 6 k 4 + 1458 N 4 k 6 r 6 - 1728 N 2 r 4 k 2 
- 259200 k 4 r 4 A 8 + 80000 k 2 r 4 A 8 + 358800 k 2 r 2 A 6 + 270000 k 2 r 4 N 4 
- 28800 r 2 A 2 - 13200 N 4 r 4 + 20000 N 6 r 4 + 99600 A 4 r 2 + 1250 N 8 r 8 k 2 
+ 56862 A 8 fc 8 r 8 + 115668 A 6 r 8 k e - 236520 A 8 fc 6 r 6 - 80550 k 4 r 8 A 8 + 720 A 2 r 4 
)/[(9 fc 2 r 2 A 2 - 3 + 5 Ar 2 )(5103 A 6 k 6 r 6 - 50625 A 6 r 6 fc 4 + 40500 A 6 k 4 r 4 

+ 10692 A 4 k 4 r 4 - 2475 A 6 k 2 r 6 - 66600 fc 2 r 4 A 6 + 37800 fc 2 r 4 N 4 
+ 50400 fc 2 r 2 A 6 - 32400 A 4 k 2 r 2 + 1296 fc 2 r 2 N 2 + 125 r 6 + 500 r 4 
- 300 N 4 r 4 - 12000 N 6 r 2 + 13200 A 4 r 2 - 3600 r 2 A 2 + 32000 A 6 - 52800 A 4 

+ 28800 A^ 2 — 5184)] 



(A.2) cos(2</>) = 3(13851 N 6 k 6 r 6 - 144585 A 6 r 6 k 4 - 10575 N 6 k 2 r 6 + 125 N 6 r 6 

+ 1 16640 N 6 k 4 r 4 - 201600 A 6 r 4 k 2 - 4000 N 6 r 4 + 176400 fc 2 r 2 - 58000 r 2 
+ 80000 A 6 + 11664 A 4 k 4 r 4 + 127440 k 2 r 4 A 4 + 2400 A 4 r 4 - 88560 k 2 r 2 N 4 
+ 73200 A 4 r 2 - 100800 A 4 - 10368 k 2 r 2 N 2 - 23040 N 2 r 2 + 34560 A 2 - 1728)r 
fe/(-1728 + 34560 N 2 - 326160 /c 2 r 2 A 4 + 864 r 2 - 7776 /c 2 r 2 + 41472 k 2 r 2 N 2 
- 90000 N 6 r 2 - 13200 N 4 r 4 - 8370 A 4 k 2 r 6 + 99600 A 4 r 2 - 100800 N 4 
- 28800 A 2 r 2 + 464400 k 2 r 2 N 6 + 20000 N 6 r 4 - 625 A 6 r 6 + 11664 fc 4 r 4 N 2 
- 7776 N 4 k 4 r 4 + 30618 N 4 k 6 r 6 + 80000 A 6 + 505440 A 6 k 4 r 4 + 152361 N 6 k 6 r 6 
+ 56475 N 6 k 2 r 6 - 702675 N 6 r 6 k 4 + 59778 A 4 r 6 A: 4 - 7776 A 2 r 4 /c 2 
+ 270000 k 2 r 4 A 4 - 460800 A 6 r 4 k 2 - 132678 N 6 r 8 k 6 + 13122 A 6 k 8 r 8 

- 450 A 6 r 8 fc 2 + 18630 A 6 r 8 k 4 + 720 A 2 r 4 + 150 N 4 r 6 ) 



LATE- TIME BEHAVIOUR OF THE TILTED BIANCHI TYPE VI h MODELS 



25 



Polynomial: 

(A.3) P(k, N, A) = -9(-24 + 130 N 2 - 125 N 4 + 25N 4 k + 135 N 4 k 3 + 54 N 4 k 4 + 54 k 2 N 2 
+ 135 N 4 k 2 - 120 N 2 fc)(-24 + 130 iV 2 - 125 TV 4 - 25 N 4 k - 135 N 4 k 3 + 54 N 4 k 4 
+ 54 k 2 N 2 + 135 iV 4 k 2 + 120 iV 2 fc)(12 + 9 k 2 N 2 - 25 N 2 ) 2 + 3 iV 2 (414720 
+ 2661 1200 N 2 + 229770000 N 6 - 132408000 A 4 - 746496 k 2 

- 403734375 A 10 k 2 - 47202750 k 6 A 10 + 6298560 A 4 k e + 40234375 A 10 
+ 550314000 k 4 N 8 - 165500000 N 8 - 38283435 N 10 k 10 + 93658275 A 10 fc 8 

- 341388000 k 2 A 6 - 28868400 N 8 k 8 - 331840800 N 8 k 6 + 115998750 N 10 k 4 
+ 725805000 N 8 k 2 + 59486400 A 4 k 4 - 7776000 k 2 N 2 + 39528000 N 4 k 2 

- 229780800 A 6 k 4 + 78732000 N 6 k 6 - 2799360 N 2 k 4 + 16533720 N 8 fc 10 
+ 4251528 N w k 12 + 20470320 A 6 k 8 )\ 2 - 10 A 4 (77760 + 3888000 N 2 

+ 10462500 A 6 - 12352500 A^ 4 — 1399680 k 2 + 9211644 TV 4 fc 8 + 25194240 A^ 4 k 6 
- 56457000 k 4 N 8 - 419904 k 4 - 1484375 A 8 + 339174000 k 2 N 6 
+ 104090265 N 8 k 8 ~ 34736850 A 8 k 6 ~ 206015625 A" 8 k 2 + 1889568 A 2 k 6 

- 61 148520 A" 4 k 4 - 2060640 k 2 N 2 - 130442400 A 4 k 2 + 414752400 N 6 k 4 

- 211789080 N 6 k 6 + 14696640 A 2 k 4 + 3188646 A" 8 fc 12 - 35724645 A 8 k w 

- 32673780 N 6 k 8 + 9920232 A 6 k w )X 4 + 30 AT 6 (7200 + 64044000 N 4 k 2 

- 19306350 N 6 k 6 - 18994095 A 6 k w - 72171000 N 4 k 6 + 1400000 AT 4 

- 2799360 A" 2 k 4 + 64921095 A 6 k 8 - 50521875 k 2 N 6 - 64739250 N 6 k 4 

- 339000 N 2 + 3306744 A 4 A: 10 - 16008840 A" 4 fc 8 + 1417176 A 6 k 12 - 375840 k 2 
+ 2047032 A" 2 fc 8 + 3324240 N 2 k 6 + 132013800 A 4 k 4 - 22366800 k 2 N 2 
+ 1283040 k 4 + 209952 k 6 - 1328125 N 6 )\ 6 - 15 Ar 8 (13599000 k 2 N 2 

- 33934950 N 4 fc 6 - 33165855 N 4 k w - 71092080 N 2 k 6 - 262500 A 2 

+ 4694760 k 4 + 130793535 A 4 k 8 - 27496875 A" 4 k 2 - 99636750 N 4 k 4 + 1500 
+ 3306744 N 2 fc 10 - 21126420 A 2 fc 8 + 2125764 A" 4 fc 12 + 1023516 fc 8 + 524880 fc 6 
+ 88176600 N 2 k 4 - 183600 k 2 + 546875 A 4 )A 8 + A" 10 (9 k 2 - 30 k + 5) 
(9 k 2 + 30 k + 5) (157464 N 2 k 8 + 122472 k 6 - 1228365 N 2 k 6 + 232875 N 2 k 4 
+ 252720 k 4 - 5400 k 2 + 410625 k 2 N 2 + 15625 A 2 )A 10 - 

324 A 12 k 4 (9 k 2 - 30 k + 5) 2 (9 k 2 + 30 k + 5) 2 A 12 

A. 2. Expressions determining and V for IZ(h). 

(A.4) fl = -(24 - 92 7 + 120 k + 168 k 2 + 72 fc 3 + 114 7 2 - 45 7 3 + 96 S+ 2 - 96 S+ - 376 k 7 
+ 366 k 2 7 2 + 366 k 7 2 - 444 k 2 7 - 93 k 2 7 3 - 108 k 7 3 - 90 k 3 7 3 + 258 k 3 j 2 
- 240 k 3 7 - 672 k 2 S + - 480 k S+ + 480 k S+ 2 + 672 k 2 S+ 2 - 288 k 3 S+ 
+ 288 fc 3 S+ 2 + 224 7 S+ - 80 S+ 2 7 - 120 7 2 S+ + 1632 k 2 S+ 7 + 276 £+ k 2 7 3 
+ 1168fc7S+ + 1200 S + 2 k 2 7 2 - 1224 £+ k 2 j 2 - 948 k 7 2 S+ - 832fc 7 S+ 2 

- 1488 k 2 7 S+ 2 - 360 S+ 2 k 2 7 3 + 270 k 7 3 S + + 360 k 7 2 S+ 2 + 126 £+ fc 3 7 3 

- 444 £+ fc 3 7 2 + 216 S+ 2 fc 3 7 3 - 216 S+ 2 fc 3 7 2 + 624 k 3 S + 7 - 288 fc 3 S+ 2 7) j [6 

7 2 k 2 (6 + 3fc 7 -5 7 -2fc)] 
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(A.5) V 2 = (24- 927+ 144fc + 216fc 2 + 114 7 2 - 457 3 + 96E + 2 -96£ + -432/c7 + 378/c 2 7 2 
+ 396 fc 7 2 - 492 fc 2 7 - 93 k 2 7 3 - 108 fc 7 3 - 90 fc 3 7 3 + 168 fc 3 7 2 - 72 fc 3 7 
- 864 k 2 £+ - 576 k £+ + 576 k S+ 2 + 864 k 2 S+ 2 + 224 7 £+ - 80 S+ 2 7 
- 120 7 2 £+ + 1968 k 2 £+ 7 + 276 £+ A: 2 7 3 + 1320 fc 7 £+ + 1440 S + 2 k 2 7 2 
- 1368 £+ fc 2 7 2 - 1008fc 7 2 S + - 912/c 7 S + 2 - 1968fc 2 7£ + 2 - 360S + 2 /c 2 7 3 
+ 270 fc 7 3 £+ + 360 fc 7 2 S+ 2 + 126 £+ fc 3 7 3 - 192 S + fc 3 7 2 + 216 S+ 2 fc 3 7 3 
- 360 S+ 2 fc 3 7 2 + 72 fc 3 £+ 7 + 144 fc 3 S + 2 7) j [ 

(6 fc 2 £+ 7 + 3 7 - 2 + 4 £+ + 5 fc 7 - 10 k 7 £+ - 6 k + 12 k S+)(60 fc 7 2 £+ 
-21fc7 2 - 157 2 + 60fc7 + 287- 132fc 7 £ + - 207X: + + 72fc£ + -36fc-12 

+ 24S+)] 
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